Libraries:

# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
library(ggeffects)
library(performance)

Set working directory

knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')

Source sapling data:

source("Scripts/SA_Import.R")

Pivot wider to create dataframe where each row is for one plot and has total details for each species group

sa_merge2 <- sapling_merge %>% 
  pivot_wider(names_from = Species_Groups, values_from = Total_Tally)

Import time since data and add it to the PIRI dataset

time_since <- read_csv("CleanData/Treat_Year_Data.csv")

sa_merge3 <- merge(sa_merge2, time_since, by = 'Site')
#log transform time from treatment data
sa_merge3$l.TFT <- log(sa_merge3$Time_from_Treat)

Run the ‘Add_BA’ script and merge with dataset:

source("Scripts/Add_BA.R")

# merge with sapling dataset -------------------
sa_merge4 <- merge(sa_merge3, prism_BA, by = 'Plot_No')

Run ‘Ground_Data.R’ script and add it to sapling dataset:

source("Scripts/Ground_Data.R")

# merge with sapling dataset -------------------
sa_merge5 <- merge(sa_merge4, ground3, by = 'Plot_No')

Source and add veg data:

source("Scripts/Veg_Data.R")

# merge with sapling dataset -------------------
sa.all <- merge(sa_merge5, veg3, by = "Plot_No")

rm(sa_merge3,
   sa_merge4,
   sa_merge5)

Sapling count data is taken in 25m2 plots; basal area is measured in hectares; veg and soil data is taken a 1m2 plots.

Sapling data will be convered into 1m2 plots in order to compare across and reduce the amount of scales of data collection, reducing it to two: 1m2 plots and per hectare observations

Create log transformed categories for newly added variables, then select for just the desired variables:

sa.all$PIRI.1m <- sa.all$PIRI/25
sa.all$SO.1m <- sa.all$Shrub_Oak/25
sa.all$Other.1m <- sa.all$Other/25

sa.all$l.PIRI1 <- log(sa.all$PIRI.1m + 1)
sa.all$l.SO1 <- log(sa.all$SO.1m + 1)
sa.all$l.other1 <- log(sa.all$Other.1m + 1)

sa.all2 <- sa.all %>% 
  select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI1, Shrub_Oak, l.SO1, Other, l.other1, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>% 
  arrange(Treat_Type)

Create a dataframe keeping sapling at 25m2 observation levels & log transforming them

sa.all3 <- sa.all

sa.all3$l.PIRI <- log(sa.all3$PIRI + 1)
sa.all3$l.SO <- log(sa.all3$Shrub_Oak + 1)
sa.all3$l.other <- log(sa.all3$Other + 1)

sa.all3 <- sa.all3 %>% 
  select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI, Shrub_Oak, l.SO, Other, l.other, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>% 
  arrange(Treat_Type)

Select just for numerical vs log and then look at paired plots for sapling transformed to 1m2 observation level:

#not log transformed (1m2)
sa.num <- sa.all2 %>% 
  select(PIRI.1m, S0.1m, Other.1m, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)

ggpairs(sa.num)
ggpairs(sa.num, aes(color = Treat_Type))


#log transformed (1m2)
sa.numl <- sa.all2 %>% 
  select(l.PIRI1, l.SO1, l.other1, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)

ggpairs(sa.numl)
ggpairs(sa.numl, aes(color = Treat_Type))

rm(sa.num,
   sa.numl)

Can see the correlation coefficients for linear (Pearsons) relationships. None of them appear very strong, except for ones that are analogs (avg LD vs mineral soil exposure; ba/ha vs piri ba/ha)

No significant correlations obvious from plots

Above code and naming conventions are the same as SA_GLMM.Rmd script

Looking at pairwise again but with sa.all3 dataset (sapling counts at 25m2 observation level)

#not log transformed (25m2)
sa3.num <- sa.all3 %>% 
  select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)

ggpairs(sa3.num)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggpairs(sa3.num, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#log transformed (25m2)
sa3.numl <- sa.all3 %>% 
  select(l.PIRI, l.SO, l.other, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)

ggpairs(sa3.numl)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggpairs(sa3.numl, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

rm(sa3.num,
   sa3.numl)

Again no strong correlations

Going to try some scatterplots:

#PIRI & SO
ggplot(sa.all3, aes(x=l.PIRI, y = l.SO))+
  geom_point()+
  geom_smooth(method = lm)+
  facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'

# negative for control, positive for fallrx, flat for other

#PIRI and other
ggplot(sa.all3, aes(x=l.PIRI, y = l.other, fill = Treat_Type))+
  geom_point()+
  geom_smooth(method = lm) +
  facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'

# negative for control & fall rx, weak negative for springrx and mow, positive for harvest

#PIRI and veg
ggplot(sa.all3, aes(x=l.PIRI, y = l.Veg_Total, fill = Treat_Type))+
  geom_point()+
  geom_smooth(method = lm) +
  facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'

# positive for harvest, weak negative for springrx, negative for mowrx, fallrx, and control

#relationships varying by treatment type may be a reason why the model are having a hard time in DHARMa

major personal breakthrough for dummies. For zero inflation, lets see if any TT or region is lacking alltogether in PIRI or SO saplings

sa.all4 <- sa.all3 %>% 
  group_by(Region)

tapply(sa.all4$PIRI, sa.all4$Region, summary)
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.02459 0.00000 2.00000 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2013  0.0000  9.0000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.07273 0.00000 1.00000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.09574 0.00000 3.00000
# no PIRI in LI,
tapply(sa.all4$PIRI, sa.all4$Treat_Type, summary) 
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2564  0.0000  9.0000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.04902 0.00000 2.00000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.05357 0.00000 2.00000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.05882 0.00000 2.00000 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01149 0.00000 1.00000

Now that I know that there are no PIRI in LI, I want to start model elimination again using ziformula = ~ Region

I have previously done model elimination with zero inflated poisson distro & with model w.o TT. Neither are better than zi nbinom2 models. zi ~1 has distribution issues, ~Region does not

Model has much better fit and zero inflation issues are resolved when treatment type variable is included. Time to do elimination again but with TT; all models fail with nbinom1 distro, but works with nbinom2 (these notes are from when I was using a zi poisson distro, using nbinom TT no longer seems to be important to model fit, but we will keep for the overarching story)

sa.m1b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
 #won't converge

# test piri ba
sa.m2 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m2) #248

check_collinearity(sa.m2) #not rank deficient

#test overall BA
sa.m2b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m2b) #239.3 (AIC smaller when PIRI BA included instead of overall BA)

check_collinearity(sa.m2b) # this model is not rank deficient

# no ba to compare with ba/piri ba models
sa.m3 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m3) #252.2

check_collinearity(sa.m3) #not rank deficent

lrtest(sa.m2b, sa.m3)  #p=0.0001, keep piri ba
lrtest(sa.m2, sa.m3) # p = 0.01

#I think these two variables will have issues of collinearity, but I would like to confirm that, but models with both won't converge or are rank deficient

#test SO
sa.m4 <- glmmTMB(PIRI ~ Treat_Type + l.other + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m4) #240

check_collinearity(sa.m4) #rank deficient
emmeans(sa.m4, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
# ok, so rank deficiency shows up in pairwise as NaN

lrtest(sa.m2b, sa.m4) # p = 0.08 # drop

#test other
sa.m5 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m5) #240.6

check_collinearity(sa.m5) #not rank deficient 

#test mineral 
sa.m6 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m6) #238.8

check_collinearity(sa.m6) #not rank deficient 

lrtest(sa.m5, sa.m6) # both other and mineral non-significant

#test veg
sa.m7 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri  + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m7) #237.8

check_collinearity(sa.m7) #RANK DEFICIENT
lrtest(sa.m6, sa.m7) #drop

sa.m7b <- glmmTMB(PIRI ~ Treat_Type + l.BA_HA  + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m7b) #247.7

check_collinearity(sa.m7b) # not rank deficient

#why is this happening
#if both BA/HA and piri BA/HA are important, should I make a proportion category? no - model using proportion won't converge


# sa.all3$BA_prop <- sa.all3$PIRI.BA_HA/sa.all3$BA_HA
# sa.all3$l.BA_prop <- log(sa.all3$BA_prop+1)
# 
# hist(sa.all3$BA_prop)
# 
# sa.m7p <- glmmTMB(PIRI ~ Treat_Type + l.BA_prop + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
#                  data = sa.all3,
#                  ziformula = ~Region,
#                  family = nbinom2)

#ok this idea failed. can't get models to work with piri ba proportion or log transformed piri ba prop

# i don't know how to resolve this issue or what the right steps to take


sa.m7c <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri  + avgLD_l + l.SO + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m7c) #237.3
# won't converge with l.other; above model is rank deficient

rm(sa.m7b, sa.m7c)







#test avgld
sa.m8 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m8) #won't converge

check_collinearity(sa.m8) #RANK DEFICIENT

#test TT
sa.m9 <- glmmTMB(PIRI ~ l.BA_piri + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m9) #238.2
lrtest(sa.m9, sa.m2b) # p = 0.06


#test piri ba
sa.m10 <- glmmTMB(PIRI ~ l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(sa.m10) #245.2



#will avgld and piri ba converge without TT? NO



# tried with zi ~1, major dispersion issues when checking model fit

rm(sa.m1b, sa.m2, sa.m2b, sa.m3, sa.m4, sa.m5, sa.m10, sa.m11)

Checking fixed effects for collinearity and looking at general summary stats to make sure they aren’t 0/1 or same/function of other category

summary(sa.all3)
##   Treat_Type           Region              Site              Plot_No      
##  Length:498         Length:498         Length:498         Min.   :  70.0  
##  Class :character   Class :character   Class :character   1st Qu.: 359.5  
##  Mode  :character   Mode  :character   Mode  :character   Median : 614.0  
##                                                           Mean   : 613.5  
##                                                           3rd Qu.: 887.5  
##                                                           Max.   :1143.0  
##       PIRI             l.PIRI          Shrub_Oak           l.SO       
##  Min.   :0.00000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :0.00000   Median :0.00000   Median :0.0000   Median :0.0000  
##  Mean   :0.09438   Mean   :0.04627   Mean   :0.2691   Mean   :0.1203  
##  3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.0000  
##  Max.   :9.00000   Max.   :2.30259   Max.   :7.0000   Max.   :2.0794  
##      Other           l.other      Time_from_Treat     l.TFT      
##  Min.   :0.0000   Min.   :0.000   Min.   : 1.0    Min.   :0.000  
##  1st Qu.:0.0000   1st Qu.:0.000   1st Qu.: 3.0    1st Qu.:1.099  
##  Median :0.0000   Median :0.000   Median : 5.0    Median :1.609  
##  Mean   :0.3514   Mean   :0.173   Mean   :10.9    Mean   :1.804  
##  3rd Qu.:0.0000   3rd Qu.:0.000   3rd Qu.: 8.0    3rd Qu.:2.079  
##  Max.   :9.0000   Max.   :2.303   Max.   :33.0    Max.   :3.497  
##      BA_HA          l.BA_HA        PIRI.BA_HA     l.BA_piri    
##  Min.   : 0.00   Min.   :0.000   Min.   : 0.0   Min.   :0.000  
##  1st Qu.: 7.00   1st Qu.:2.079   1st Qu.: 5.0   1st Qu.:1.792  
##  Median :16.00   Median :2.833   Median :11.0   Median :2.485  
##  Mean   :16.47   Mean   :2.521   Mean   :12.6   Mean   :2.218  
##  3rd Qu.:23.00   3rd Qu.:3.178   3rd Qu.:18.0   3rd Qu.:2.944  
##  Max.   :51.00   Max.   :3.951   Max.   :44.0   Max.   :3.807  
##   Mineral_Soil       l.Mineral       Litter_Duff        avgLD       
##  Min.   : 0.0000   Min.   :0.0000   Min.   : 0.00   Min.   : 0.000  
##  1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.:90.50   1st Qu.: 2.250  
##  Median : 0.0000   Median :0.0000   Median :90.50   Median : 3.875  
##  Mean   : 0.7199   Mean   :0.2032   Mean   :88.78   Mean   : 4.179  
##  3rd Qu.: 0.0000   3rd Qu.:0.0000   3rd Qu.:90.50   3rd Qu.: 5.500  
##  Max.   :50.5000   Max.   :3.9416   Max.   :90.50   Max.   :13.250  
##     avgLD_l        Veg_Total       l.Veg_Total   
##  Min.   :0.000   Min.   :  4.00   Min.   :1.386  
##  1st Qu.:1.179   1st Qu.: 37.50   1st Qu.:3.624  
##  Median :1.584   Median : 58.00   Median :4.060  
##  Mean   :1.519   Mean   : 62.64   Mean   :3.970  
##  3rd Qu.:1.872   3rd Qu.: 83.38   3rd Qu.:4.423  
##  Max.   :2.657   Max.   :202.00   Max.   :5.308
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
sa_check <- sa.all3 %>% 
  select(l.PIRI, l.other, l.TFT, l.BA_piri, l.Mineral, l.Veg_Total, avgLD_l)

findLinearCombos(sa.all3[,4:23]) #shows no linear combos. Don't get how the rank deficiency is happening
## $linearCombos
## list()
## 
## $remove
## NULL
tapply(sa.all3$Other, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.3443  0.0000  6.0000 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.411   0.000   6.000 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2143  0.0000  6.0000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.6727  0.0000  9.0000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.3511  0.0000  4.0000
tapply(sa.all3$Shrub_Oak, sa.all3$Region, summary) #none in LI
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $MA
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   1.964   3.000   7.000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2447  0.0000  6.0000
tapply(sa.all3$BA_HA, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   11.00   14.98   23.00   51.00 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   11.00   21.00   19.81   28.00   39.00 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    9.00   14.00   15.32   21.00   44.00 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   13.50   23.00   22.13   31.00   44.00 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   11.00   14.37   21.00   48.00
tapply(sa.all3$PIRI.BA_HA, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       2       9      10      16      41 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   11.00   16.00   17.41   23.00   37.00 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00    9.00   11.44   16.00   34.00 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    8.00   18.00   19.18   30.00   44.00 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00    9.00   10.32   16.00   32.00
tapply(sa.all3$Veg_Total, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.00   51.50   71.00   70.96   91.88  166.50 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.00   23.00   39.00   44.52   61.00  128.00 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    6.50   33.75   52.75   56.93   74.38  161.50 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   14.50   41.25   64.00   62.97   79.75  114.00 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.00   46.00   62.25   75.09   96.62  202.00
tapply(sa.all3$Mineral_Soil, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.586   0.500  30.500 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.589   0.500  15.500 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.6623  0.0000 50.5000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01818 0.00000 0.50000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2021  0.0000  8.0000
tapply(sa.all3$avgLD, sa.all3$Region, summary) #all present
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   1.750   1.977   2.750   7.250 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.75    4.00    5.25    5.61    6.75   12.50 
## 
## $MA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.250   3.562   5.000   5.347   6.750  13.250 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.750   2.500   3.500   3.686   4.500   9.750 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.750   2.812   4.000   4.298   5.438  12.250
tapply(sa.all3$Other, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.8889  1.0000  6.0000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.4118  0.0000  9.0000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1429  0.0000  2.0000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.06618 0.00000 2.00000 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1379  0.0000  2.0000
tapply(sa.all3$Shrub_Oak, sa.all3$Treat_Type, summary) #none in mowrx or springrx
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.7094  0.0000  7.0000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.4902  0.0000  6.0000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0
tapply(sa.all3$BA_HA, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   18.00   28.00   26.78   37.00   51.00 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    9.00   16.00   17.44   23.00   44.00 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   11.00   12.11   18.00   34.00 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   7.000   9.059  14.000  37.000 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   10.00   14.00   15.85   21.00   41.00
tapply(sa.all3$PIRI.BA_HA, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    9.00   18.00   17.05   25.00   41.00 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   11.00   14.45   20.25   44.00 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   5.000   9.000   9.804  14.500  30.000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   7.000   7.816  11.000  37.000 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    9.00   14.00   13.75   18.00   34.00
tapply(sa.all3$Veg_Total, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.00   26.50   46.00   46.85   61.00  115.50 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    14.5    46.0    68.0    72.7    91.0   202.0 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.00   36.12   55.00   55.21   68.62  132.00 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.00   37.00   66.75   70.47   94.38  196.50 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.00   49.50   61.50   64.65   81.00  131.00
tapply(sa.all3$Mineral_Soil, sa.all3$Treat_Type, summary) #all present
## $Control
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.000000 0.000000 0.008547 0.000000 0.500000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.201   0.000   8.000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.652   0.000  50.500 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.199   0.500  30.500 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.9368  0.5000 15.5000
tapply(sa.all3$avgLD, sa.all3$Treat_Type, summary) #all present
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.500   3.250   4.750   5.024   6.750  12.500 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.750   3.250   4.375   4.529   5.750   9.750 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.250   3.688   4.775   5.291   6.750  11.000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.500   2.750   2.816   3.812  13.250 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.500   1.500   3.500   4.046   5.250  11.500

Best model (lowest AIC without rank deficiency) appears to be model sa.m6 with TT, piri BA, agvLD, Veg total. Test model fit:

sa.m6 <- glmmTMB(PIRI ~ Treat_Type + l.BA_piri + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
summary(sa.m6) #238.8
##  Family: nbinom2  ( log )
## Formula:          
## PIRI ~ Treat_Type + l.BA_piri + avgLD_l + l.Veg_Total + offset(l.TFT) +  
##     (1 | Site/Plot_No)
## Zero inflation:        ~Region
## Data: sa.all3
## 
##      AIC      BIC   logLik deviance df.resid 
##    238.8    306.2   -103.4    206.8      482 
## 
## Random effects:
## 
## Conditional model:
##  Groups       Name        Variance      Std.Dev. 
##  Plot_No:Site (Intercept) 0.00000002537 0.0001593
##  Site         (Intercept) 0.58921621688 0.7676042
## Number of obs: 498, groups:  Plot_No:Site, 498; Site, 47
## 
## Dispersion parameter for nbinom2 family (): 0.551 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -5.4666     2.1272  -2.570  0.01017 *  
## Treat_TypeFallRx    -0.1423     0.8538  -0.167  0.86759    
## Treat_TypeHarvest    1.6832     1.1041   1.524  0.12740    
## Treat_TypeMowRx      1.6446     0.8920   1.844  0.06522 .  
## Treat_TypeSpringRx  -0.5944     1.3264  -0.448  0.65403    
## l.BA_piri            1.3397     0.3983   3.363  0.00077 ***
## avgLD_l             -0.7821     0.6039  -1.295  0.19533    
## l.Veg_Total         -0.3636     0.3773  -0.964  0.33516    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     1.9116     0.9406   2.033   0.0421 *
## RegionLI       21.4369 30910.2106   0.001   0.9994  
## RegionMA      -20.0122 13973.6041  -0.001   0.9989  
## RegionME       -1.7758     1.3951  -1.273   0.2031  
## RegionNH       -1.5524     1.1818  -1.314   0.1890  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sa.m6_sr <- simulateResiduals(sa.m6, n = 1000, plot = TRUE) #looks good

testResiduals(sa.m6_sr) #passes

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.046286, p-value = 0.2364
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.85775, p-value = 0.688
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.046286, p-value = 0.2364
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.85775, p-value = 0.688
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 ) 
##                                                  0
testZeroInflation(sa.m6_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99994, p-value = 1
## alternative hypothesis: two.sided

Pairwise:

emmeans(sa.m6, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
##  Treat_Type response     SE  df asymp.LCL asymp.UCL
##  Control      0.0361 0.0286 Inf   0.00765     0.170
##  FallRx       0.0313 0.0253 Inf   0.00641     0.153
##  Harvest      0.1942 0.1760 Inf   0.03288     1.147
##  MowRx        0.1869 0.1217 Inf   0.05215     0.670
##  SpringRx     0.0199 0.0247 Inf   0.00175     0.227
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast           ratio     SE  df null z.ratio p.value
##  Control / FallRx   1.153  0.984 Inf    1   0.167  0.9998
##  Control / Harvest  0.186  0.205 Inf    1  -1.524  0.5465
##  Control / MowRx    0.193  0.172 Inf    1  -1.844  0.3483
##  Control / SpringRx 1.812  2.403 Inf    1   0.448  0.9917
##  FallRx / Harvest   0.161  0.186 Inf    1  -1.579  0.5106
##  FallRx / MowRx     0.167  0.163 Inf    1  -1.839  0.3510
##  FallRx / SpringRx  1.572  2.132 Inf    1   0.333  0.9973
##  Harvest / MowRx    1.039  1.138 Inf    1   0.035  1.0000
##  Harvest / SpringRx 9.753 14.468 Inf    1   1.535  0.5393
##  MowRx / SpringRx   9.384 12.813 Inf    1   1.640  0.4718
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log scale

None significantly different from eachother

Plot model

sa.piri1 <- ggpredict(sa.m6, terms = c("avgLD_l", "Treat_Type"))
#control and fallrx are almost completely identical - if I wanted to nudge/jitter
# sa.piri2 <- sa.piri1 %>% 
#   filter(group == "FallRx")
# 
# sa.piri2$predicted <- sa.piri2$predicted+0.01
# 
# sa.piri3 <- sa.piri1 %>% 
#   filter(group != "FallRx")
# 
# sa.piri3 <- rbind(sa.piri3, sa.piri2)

ggplot(sa.piri1, aes(x, predicted, color = group))+
  geom_line(linewidth = 1, show.legend = FALSE) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(axis.text = element_text(size = 10), axis.title = element_text(size = 14))+
  labs(x = "Average leaf litter depth (log transformed)",
       y = NULL)

ggsave(filename = "SA_PIRI_avgld.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)






sa.piri4 <- ggpredict(sa.m6, terms = c("l.BA_piri", "Treat_Type"))
#again control and fall rx overlap - if I wanted to nudge/jitter

# sa.piri5 <- sa.piri4 %>% 
#   filter(group == "FallRx")
# 
# sa.piri5$predicted <- sa.piri5$predicted+0.01
# 
# sa.piri6 <- sa.piri4 %>% 
#   filter(group != "FallRx")
# 
# sa.piri6 <- rbind(sa.piri6, sa.piri5)


ggplot(sa.piri4, aes(x, predicted, color = group))+
  geom_line(linewidth = 1, show.legend = FALSE) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(axis.text = element_text(size = 10), axis.title = element_text(size = 14))+
  labs(x = "Average pitch pine basal area per hectare (log transformed)",
       y = "Pitch pine saplings \n (corrected for time)")

ggsave(filename = "SA_PIRI_BA.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)





sa.piri5 <- ggpredict(sa.m6, terms = c("l.Veg_Total", "Treat_Type"))

ggplot(sa.piri5, aes(x, predicted, color = group))+
  geom_line(linewidth = 1, show.legend = FALSE) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(axis.text = element_text(size = 10), axis.title = element_text(size = 14))+
  labs(x = "Average understory vegetation cover (log transformed)",
       y = NULL)

ggsave(filename = "SA_PIRI_veg.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)

Trying to plot with jtools package

Shrub Oak Models

Start from the beginning knowing that both Region and Treat Type have 0 areas

tapply(sa.all4$Shrub_Oak, sa.all4$Region, summary)
## $ALB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000 
## 
## $LI
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $MA
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000 
## 
## $ME
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   1.964   3.000   7.000 
## 
## $NH
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2447  0.0000  6.0000
# no SO in LI,
tapply(sa.all4$Shrub_Oak, sa.all4$Treat_Type, summary) 
## $Control
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.7094  0.0000  7.0000 
## 
## $FallRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.4902  0.0000  6.0000 
## 
## $Harvest
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000 
## 
## $MowRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0 
## 
## $SpringRx
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0
# no SO in MowRX or SpringRx

rm(sa.all4)

testing with nbinom1 distro

so.sa5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')
#won't converge

so.sa6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.sa6) #334.7 - this is the only way the model will converge, to remove Mineral
## [1] 334.7261
lrtest(so.sa5, so.sa6)
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + 
##     l.BA_piri + l.BA_HA + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.PIRI + l.other + avgLD_l + l.BA_piri + 
##     l.BA_HA + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1  19                            
## 2  18 -149.36 -1
so.sa7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + avgLD_l + l.BA_piri + l.BA_HA +  offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.sa7) #333.3
## [1] 333.3224
so.sa8 <- glmmTMB(Shrub_Oak ~ Treat_Type + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.sa8) #331.3
## [1] 331.3264
so.sa9 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + avgLD_l + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.sa9) #326.3 (finally able to add l.Mineral back in)
## [1] 326.2949
lrtest(so.sa8, so.sa9) #mineral very significant
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.Mineral + avgLD_l + l.BA_HA + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq            Pr(>Chisq)    
## 1  16 -149.66                                    
## 2  16 -147.15  0 5.0315 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.sa10 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
AIC(so.sa10) #324.3
## [1] 324.3074
lrtest(so.sa9, so.sa10)
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + l.Mineral + avgLD_l + l.BA_HA + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.Mineral + l.BA_HA + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  16 -147.15                     
## 2  15 -147.15 -1 0.0125     0.9111
so.sa11 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom1)
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')
AIC(so.sa11) #won't converge
## [1] 324.2177

Test Model Fit

so.sa10_sr <- simulateResiduals(so.sa10, n = 1000, plot = T) #passes

testResiduals(so.sa10_sr) #passes

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.034621, p-value = 0.5893
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.87777, p-value = 0.77
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002008032
## sample estimates:
## outlier frequency (expected: 0.000542168674698795 ) 
##                                                   0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.034621, p-value = 0.5893
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.87777, p-value = 0.77
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002008032
## sample estimates:
## outlier frequency (expected: 0.000542168674698795 ) 
##                                                   0
testZeroInflation(so.sa10_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99318, p-value = 0.654
## alternative hypothesis: two.sided

Pairwise

emmeans(so.sa10, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
##  Treat_Type  response    SE  df            asymp.LCL asymp.UCL
##  Control    0.0000155 0.104 Inf 0.000000000000000222       Inf
##  FallRx     0.0000390 0.263 Inf 0.000000000000000222       Inf
##  Harvest    0.0000818 0.552 Inf 0.000000000000000222       Inf
##  MowRx      0.0000000 0.000 Inf 0.000000000000000222       Inf
##  SpringRx   0.0000000 0.000 Inf 0.000000000000000222       Inf
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                 ratio              SE  df null z.ratio p.value
##  Control / FallRx             0               0 Inf    1  -2.621  0.0665
##  Control / Harvest            0               0 Inf    1  -1.415  0.6178
##  Control / MowRx     2699054976  79659759310326 Inf    1   0.001  1.0000
##  Control / SpringRx    37242947    314408982995 Inf    1   0.002  1.0000
##  FallRx / Harvest             0               1 Inf    1  -0.624  0.9713
##  FallRx / MowRx      6809895316 200986873784945 Inf    1   0.001  1.0000
##  FallRx / SpringRx     93966433    793274786998 Inf    1   0.002  1.0000
##  Harvest / MowRx    14295160591 421906579510505 Inf    1   0.001  1.0000
##  Harvest / SpringRx   197251968   1665222450358 Inf    1   0.002  1.0000
##  MowRx / SpringRx             0             424 Inf    1   0.000  1.0000
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log scale

Plot model

so.p1 <- ggpredict(so.sa10, terms = c("l.BA_HA", "Treat_Type"))

ggplot(so.p1, aes(x, predicted, color = group))+
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
  labs(x = "Average basal area per hectare (log transformed)",
       y = "Total Count of Shrub Oak",
       title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH")

#very small effect; springrx and mowrx both have 0 SO saplings





so.p2 <- ggpredict(so.sa10, terms = c("l.Mineral", "Treat_Type"))

ggplot(so.p2, aes(x, predicted, color = group))+
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
  labs(x ="Average exposed mineral soil (log transformed)",
       y = "Total Count of Shrub Oak",
       title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH") +
  xlim(0, 0.5)
## Warning: Removed 25 rows containing missing values or values outside the scale range
## (`geom_line()`).

this is a weird one, because exposed mineral soil clearly doesn’t go over a certain percent; the median is 0 and the mean 0.2, the IQR 0 - 0. Maybe it looks better at a smaller x scale.

summary(sa.all3$l.Mineral)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2032  0.0000  3.9416
summary(sa.all3$Mineral_Soil)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.7199  0.0000 50.5000
I played around a little with using Treat Type (instead of and with Region) for zi, but it didn’t work
so.sa1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + l.BA_piri + l.BA_HA + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(so.sa1) #333.8
## [1] 333.7745
# check_collinearity(so.sa1) #rank deficient

# test avgld
so.sa2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_piri + l.BA_HA + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(so.sa2) #331.9
## [1] 331.8841
check_collinearity(so.sa2) #not rank deficient
## # Check for Multicollinearity
## 
## * conditional component:
## 
## Low Correlation
## 
##        Term  VIF       VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##  Treat_Type 1.06 [1.01,     1.30]         1.03      0.95     [0.77, 0.99]
##   l.BA_piri 2.38 [2.10,     2.74]         1.54      0.42     [0.37, 0.48]
##     l.BA_HA 2.48 [2.18,     2.86]         1.58      0.40     [0.35, 0.46]
##      l.PIRI 1.06 [1.01,     1.30]         1.03      0.95     [0.77, 0.99]
##     l.other 1.23 [1.13,     1.39]         1.11      0.81     [0.72, 0.88]
##   l.Mineral 1.00 [1.00,      Inf]         1.00      1.00     [0.00, 1.00]
# test piri ba
so.sa3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_HA + l.PIRI + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(so.sa3) #329.9
## [1] 329.9233
lrtest(so.sa2, so.sa3) #drop
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + l.BA_piri + l.BA_HA + l.PIRI + l.other + 
##     l.Mineral + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.BA_HA + l.PIRI + l.other + l.Mineral + 
##     offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  18 -147.94                     
## 2  17 -147.96 -1 0.0393     0.8429
# test piri
so.sa4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.BA_HA + l.other + l.Mineral + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = nbinom2)
AIC(so.sa4) #328
## [1] 327.9643
lrtest(so.sa3, so.sa4)
## Likelihood ratio test
## 
## Model 1: Shrub_Oak ~ Treat_Type + l.BA_HA + l.PIRI + l.other + l.Mineral + 
##     offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + l.BA_HA + l.other + l.Mineral + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1  17 -147.96                    
## 2  16 -147.98 -1 0.041     0.8396
# anything dropped from model won't converge, higher aic with ~1, won't converge with ~TT
check_collinearity(so.sa4) #not rank deficient
## # Check for Multicollinearity
## 
## * conditional component:
## 
## Low Correlation
## 
##        Term  VIF       VIF 95% CI Increased SE Tolerance Tolerance 95% CI
##  Treat_Type 1.02 [1.00,     4.31]         1.01      0.98     [0.23, 1.00]
##     l.BA_HA 1.07 [1.02,     1.28]         1.04      0.93     [0.78, 0.98]
##     l.other 1.07 [1.02,     1.28]         1.03      0.94     [0.78, 0.98]
##   l.Mineral 1.00 [1.00,      Inf]         1.00      1.00     [0.00, 1.00]

Check model fit

so.sa4_sr <- simulateResiduals(so.sa4, n = 1000, plot = TRUE) #looks good

testResiduals(so.sa4_sr) #passes

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.031283, p-value = 0.7144
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.8769, p-value = 0.93
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000542168674698795 ) 
##                                                   0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.031283, p-value = 0.7144
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.8769, p-value = 0.93
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000542168674698795 ) 
##                                                   0
testZeroInflation(so.sa4_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99681, p-value = 0.902
## alternative hypothesis: two.sided

Pairwise

emmeans(so.sa4, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
##  Treat_Type   response     SE  df            asymp.LCL asymp.UCL
##  Control    0.00000260 0.0935 Inf 0.000000000000000222       Inf
##  FallRx     0.00000602 0.2162 Inf 0.000000000000000222       Inf
##  Harvest    0.00000958 0.3439 Inf 0.000000000000000222       Inf
##  MowRx      0.00000000 0.0000 Inf 0.000000000000000222       Inf
##  SpringRx   0.00000000 0.0000 Inf 0.000000000000000222       Inf
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                 ratio               SE  df null z.ratio p.value
##  Control / FallRx             0                0 Inf    1  -2.149  0.1996
##  Control / Harvest            0                0 Inf    1  -1.004  0.8539
##  Control / MowRx     7591553109  315497435464353 Inf    1   0.001  1.0000
##  Control / SpringRx    42869891     323998497992 Inf    1   0.002  1.0000
##  FallRx / Harvest             1                1 Inf    1  -0.353  0.9967
##  FallRx / MowRx     17556161776  729616711771682 Inf    1   0.001  1.0000
##  FallRx / SpringRx     99140550     749276197656 Inf    1   0.002  1.0000
##  Harvest / MowRx    27922185717 1160418410274454 Inf    1   0.001  1.0000
##  Harvest / SpringRx   157678021    1191685840742 Inf    1   0.002  1.0000
##  MowRx / SpringRx             0              239 Inf    1   0.000  1.0000
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log scale

None significantly different. No SO saplings in SpringRx or MowRx - not sure if I should deal with that differently?

now to plot

sa.so1 <- ggpredict(so.sa4, terms = c("l.BA_HA", "Treat_Type"))

# if I wanted to jitter/nudge lines
# sa.so2 <- sa.so1
# sa.so2 <- sa.so2 %>% 
#   filter(group == "SpringRx")
# 
# sa.so2$predicted <- sa.so2$predicted+0.05
# 
# sa.so3 <- sa.so1 %>% 
#   filter(group != "SpringRx")
# 
# sa.so3 <- rbind(sa.so3, sa.so2)


ggplot(sa.so1, aes(predicted,x, color = group))+
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  coord_flip()+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
  labs(y = "Average basal area per hectare (log transformed)",
       x = "Total Count of Shrub Oak",
       title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH")

# such a small effect








sa.so2 <- ggpredict(so.sa4, terms = c("l.Mineral", "Treat_Type"))

ggplot(sa.so2, aes(predicted,x, color = group))+
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  coord_flip()+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
  labs(y = "Average exposed mineral soil per 1m^2^ (log transformed)",
       x = "Total Count of Shrub Oak",
       title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH")

#VERY WEIRD, likely because mineral soil is so limited



sa.so3 <- ggpredict(so.sa4, terms = c("l.other", "Treat_Type"))

ggplot(sa.so3, aes(predicted,x, color = group))+
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#D8B70A",  "#02401B", "#A2A475", "#81A88D", "#972D15"))+
  coord_flip()+
  theme_classic()+
  theme(panel.background = element_blank()) +
  labs(colour = "Treatment Type")+
  theme(legend.position = "right", legend.text = element_text(size = 12), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size = 14), plot.title=element_text(hjust=0.5, size=20))+
  labs(y = "Total counts of other sapling species (log transformed)",
       x = "Total Count of Shrub Oak",
       title = "Predicted Counts of Shrub Oak Saplings \n >/= 2.5 cm & < 10 cm DBH")

# so small